%This program simulates the equilibrium paths leading to the Ponzi steady
%state.
%Produces Figure 3
tic
clear;
format long;

T=350; %Number of years of horizon

%% Calibration

%Exogenous parameters
rho=0.04; %Discount rate
n=-0.007*0-0.017; %Population growth rate
g=0.013; %Rate of technical progress
frisch=0.85; %Frisch elasticity of labor supply
piW=g; %Downward wage rigidity

%Targets to match
LaborSupply=1; %Labor supply in (neoclassical) steady state
ConsumptionGap=0.1; % -(cSS-cFB)/cFB
MaxPonzi=1.25; %Ratio of Ponzi scheme to GDP in the Ponzi steady state
HalfLife=2; %Years necessary for the inflation anchor to close half the gap with the (constant) rate of inflation
PhillipsSlope=-0.3; %Elasticity of nominal wage inflation with respect to unemployment
MinWealth=2; %(Negative of) minimum wealth relative to consumption under secular stagnation

%Endogenous parameters
kL=LaborSupply^(-1-1/frisch); %Scale parameter of labor supply function
theta=log(2)/HalfLife; %Speed of inflation deanchoring
beta=-PhillipsSlope*frisch; %Slope coefficient of the Phillips curve

%Neoclassical steady state & Secular stagnation steady state
LNeo=LaborSupply;
LSS=(1-ConsumptionGap)*LNeo;
NaturalRate=g-rho*(LNeo/LSS-1)-LNeo/LSS*piW;
lSS=(kL*LSS)^(-frisch); %Labor supply

%Ponzi steady state
deltaPonzi=MaxPonzi*LNeo;

%Calibration of the preference for wealth
RefWealth=-MinWealth*LSS;
sigma=log(LNeo*(rho+piW)/(LSS*(rho-n)))/log(1-deltaPonzi/RefWealth);
kW=((rho-n)/LNeo)*(deltaPonzi-RefWealth)^sigma;

%Taylor rule
pix=0.02; %Inflation target
phi=1.5; %Taylor parameter

%Public debt
S=10; %Maturity of long-term debt
D=1.66; %Long-term indebtedness
H=0.25*LSS*2; %Helicopter drop

%% Simulating the equilibrium trajectory that is plotted

%Initial values
piA0=piW-g;
delta0=H;
L0=0.91340890714; %This value was obtained by hand (this is the value of L0 such that delta tends to MaxPonzi).

%Solves a differential equation from time 0 to T using ode15s or ode23t
sol=ode15s(@(t,y) Equi(t,y,rho,g,n,frisch,kL,theta,beta,phi,pix,kW,RefWealth,sigma,NaturalRate,piW),[0 T],[L0;piA0;delta0;0]);
t=sol.x';
y=sol.y';
%yT=deval(sol,T);
y(end,3);

L=interp1(t,y(:,1),0:0.01:T)';
piA=interp1(t,y(:,2),0:0.01:T)';
delta=interp1(t,y(:,3),0:0.01:T)';
pi=max(piA+beta*(kL*L.^(1+1/frisch)-1),piW-g);
i=max(NaturalRate+pix+phi*(pi-pix),0);

t=[0:0.01:T]';

%% Figures

FontMain=34;
Font=26;
LW=3.4;
ST=-20;
TT=T;

%Consumption
axes('Parent',figure('InvertHardcopy','off','Color',[1 1 1]),'box','on','FontSize',Font,'position',[0.04 0.58 0.43 0.38],'LineWidth',1)
hold on
xlim([ST TT])
ylim([0.89 1.02])
plot([ST;0;t(1:TT*100)],[LSS;LSS;L(1:TT*100,1)],'DisplayName','Ponzi','LineWidth',LW,'LineStyle','-','Color',[0 0 0])
xlabel('Years','FontSize',Font)
title('Output','FontSize',FontMain)

%Inflation
axes('box','on','FontSize',Font,'position',[0.53 0.58 0.43 0.38],'LineWidth',1)
hold on
xlim([ST TT])
ylim([-0.1 0.55])
plot([ST;0;t(1:TT*100)],[100*(piW-g);100*(piW-g);100*piA(1:TT*100,1)],'LineWidth',LW,'LineStyle','-','Color',[0 0 0])
%plot(t(1:TT*100),100*pi(1:TT*100,1),'LineWidth',1.5,'LineStyle','-','Color',[0 0 0])
xlabel('Years','FontSize',Font)
title('Inflation anchor (%)','FontSize',FontMain)

%Ponzi
axes('box','on','FontSize',Font,'position',[0.04 0.07 0.43 0.38],'LineWidth',1)
hold on
xlim([ST TT])
ylim([-0.08 1.45])
plot([ST;0;t(1:TT*100)],[0;0;delta(1:TT*100,1)],'LineWidth',LW,'LineStyle','-','Color',[0 0 0])
xlabel('Years','FontSize',Font)
title('Ponzi scheme','FontSize',FontMain)

%Nominal interest rate
axes('box','on','FontSize',Font,'position',[0.53 0.07 0.43 0.38],'LineWidth',1)
hold on
xlim([ST TT])
ylim([-0.5 2])
plot([ST;0;t(1:TT*100)],[0;0;100*i(1:TT*100,1)],'LineWidth',LW,'LineStyle','-','Color',[0 0 0])
xlabel('Years','FontSize',Font)
title('Nominal interest rate (%)','FontSize',FontMain)

toc
